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Abstract. To simulate open boundaries within finite computation domain, real- 
function coordinate transformation in the framework of generally covariant formulation 
of Maxwell equations is proposed. The mapping — realized with arctangent function 
here — has a transparent geometric meaning of pure squeezing of space, is admissible 
by classical electrodynamics, does not introduce artificially lossy layers (or 'lossy 
coordinates') to absorb outgoing radiation nor leads to non-Maxwellian fields. At the 
same time, like for anisotropic perfectly matched layers, no modification (except for 
transformation of material tensors) is needed to existing nearest-neighbor computation 
schemes, which makes it well suited for parallel computing implementation. 



1. Introduction 

Direct numerical methods of modeling electromagnetic phenomena, such as finite- 
difference time-domain (FDTD) and frequency- domain (FDFD) schemes, are invariably 
concerned with how to represent infinite space surrounding the region of interest on a 
bounded computation domain. Two approaches to combat that problem do exist; one, 
as in |1, or f^, can be classified as non-local since it is not limited to the treatment 
of nearest-neighbor interactions on a computation grid due to higher-order differentials 
appearing in the formulation, a shortcoming when parallel computer implementation 
is considered; another is aimed to modify the (local) material properties of boundary 
regions in such a way that outgoing radiation experiences no parasitic refiections from 
the boundaries of computation window OH], and hence can be called local. 

The now-classic local technique to represent open boundaries (called absorbing 
boundaries when zero refiectivity of surrounding space is emphasized and mimicked 
in simulations) is by means of absorbing perfectly matched layers (PMLs) introduced 
originally in The technique, especially in its non-split version is considered simple 
and efficient, and enjoy great popularity in the electromagnetic modeling community. 
Nonetheless it cannot be considered as completely perfect because: 

(i) free PML parameters, such as maximum conductivity and conductivity profile, bear 
not always obvious geometric or physical relation to particular problem to simulate, 
hence need for the adjustment and optimization of PMLs so often; 
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(ii) use of complex-valued matrices for modified dielectric permittivity e and magnetic 
permeability fi leaves no room for CPU time and memory savings with real-field 
finite-difference formulations; 

(iii) loss of accuracy when solving eigenproblems in frequency domain is inevitable 
(though normally minor) due to nonzero mode tails extending into the regions 
of modified e and /i within PMLs; 

(iv) difficulties with non-Cartesian and, in particular, non-orthogonal grids have been 
reported 

As an attempt to overcome these problems while retaining locality of the formulation, 
we propose a conceptually simple and numerically easy-to-implement squeezing of open 
boundaries (SOB) technique in this Letter. 

The underlying idea of SOB is to map infinite surrounding space (or rather 
the whole space, with better sampling for central region and coarser for outskirts) 
onto the finite computation domain, instead of inserting anisotropic absorbing PMLs 
between the region of interest and computation boundaries. A clever way to do such 
mapping inexpensively is by transforming e and fi fields as stipulated by generally 
covariant electrodynamics, while retaining the form of Maxwell equations untouched. 
The mapping — illustrated by use of arctangent function here, with other possibilities 
discussed — is rigorous at the stage of analytic description; is smooth, while the nth 
derivative of material tensors at PML interface (with n depending on the order of the 
profile) is discontinuous; and is real-valued, enabling finite-difference algorithms in real 
notation where appropriate. Another advantage of the SOB method is its extendability: 
that is, anisotropic, magnetic materials and nontrivial backgrounds can be treated 
straightforwardly; and with same ease non-Cartesian and non-orthogonal coordinates, 
if preferred, can be squeezed in the manner proposed. Finally, this technique justifies a 
surprising possibility for lossless PML formulation. 

2. Covariant Maxwell equations 

It was Lorentz-covariance of Maxwell equations that led to special relativity over a 
century ago. Another, less celebrated though well and long ago established fact about 
Maxwell equations is that they can be formulated in a generally covariant manner, 
i.e., so that they do not change their form under arbitrary reversible transformation 
from Cartesian coordinates jZllH]. Surprisingly, this feature was first exploited in direct 
computation electromagnetics perhaps only a decade ago, in 'logically Cartesian' FDTD 
simulations of high index contrast dielectric structures [01 CH] (see also fij). For the 
sake of completeness let us write coordinate-invariant Maxwell equations here, in terms 
of electric covariant vector E,y and magnetic covariant pseudovector Hi, as in ^2] : 




(2) 
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a form which, written in components exphcitly, is identical to conventional Cartesian 
representation (pseudo permutation field fB'^^'^ equals Levi-Civita symbol in any 
coordinate system), with the constitutive relations 

= /'^iJ^, S)^ = e^-'E^ (3) 

(i.e., no optical activity assumed), where is magnetic induction pseudo vector 
density of weight +1, and 2)^ is electric induction vector density; hence e'^^ and fi^'^ 
are contravariant tensor densities transformed according to 

e^" = |ArV^,j;,e^''^', = |ArV^,J,n/''^', (4) 

where Jy = dyx^ is the Jacobian transformation matrix for contravariant components, 
A = det Jy is its determinant. 

Such formulation enables one to hide all metric information into e''^" and fi''^" while 
invariably using Cartesian-like representation of Maxwell equations (P), (j21), which is 
extremely convenient; no wonder many authors strived (successfully) to 'derive' (jlj) on 
different grounds and under different assumptions, as in for isotropic media or in |Tn] 
for diagonal transformation matrices Jy. It is worth noting that, however, are a 
direct consequence of transformation characteristics assigned to electric and magnetic 
fields; hence their general nature and no need in any tricky derivations. In practice, 
dielectric permittivity and magnetic permeability are referenced to Cartesian frame, so 
if one wants to use non-Cartesian coordinates instead, with Cartesian-like equations (^, 
((2)), then transformation rules (ji} are to be employed and specified for transformation 
from Cartesian to those coordinates. We make such specification for the mapping onto 
arctangent-squeezed coordinates in the next Section. 

3. Open boundaries on confined domain 

Let the Cartesian coordinates {x^' , x^\ x^'} = {x, y, z} be transformed to {x^, x^, x^} = 
{u,v,w} according to 

u = arctan(x/xo), v = arctan(?//?/o), w = aictaii{z / zq) . (5) 

where xq, yo, zq are the units of length along the corresponding coordinates. Such 
mapping preserves central region of space (where the scatterer or the waveguide is 
supposedly located) virtually untouched, while smoothly squeezing the outer space 
into the (— 7r/2, 7r/2)'^ bounded computation domain. The xq, yo and units are 
arbitrary at this, analytic stage, but their careless choice may compromise accuracy of 
finite-difference calculations on (more or less) equidistant grids in squeezed coordinates. 
Indeed, poor sampling on a 'squeezed' grid of 'physical' lengths far from the origin of 
coordinates (|x| ^ xq, \y\ ^ yo, or \A ^ -^o) is a reason to match scaling factors xq, 
and Zq with actual physical dimensions of the region of interest or with the wavelength. 
Note that in the proposed approach, xq, yo and z^ are the only parameters to be adjusted 
to particular physical problem, with a clear geometric relation to that problem. 
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By differentiating (0) one gets the transformation matrix 

/ Sg^COS^M \ 

yQ^ cos^ V 



(^0 Zq'^cos'^w j 



(6) 



The 'squeezed' permittivity and permeabihty can be obtained with (jH), © immediately; 
the diagonal components of ^ for example, are 

^uu ^ ^xx yO^O COS^ U ^ ZqXq COS^ V ^ XqI/q COS^ W 

Xo COS^ V COS^ 1(7 ' yo COS^ M COS^ W ' Zo cos^ U COS^ V ' 

and similarly for fi^^. It should be emphasized that even if matrix representations 
of dielectric permittivity and magnetic permeability are non-diagonal in Cartesian 
coordinates, this generates no additional complexity in deriving the e'^^ and fi^'^ off- 
diagonal components; here we omit them for brevity. Specifying the transformed e'*''' 
and /i'*'^ for squeezed cylindrical or spherical coordinates also poses no difficulty in 
present approach, but is out of scope here. 

The structure of (|7j) resembles that of modified e and fi in PML regions with 
'stretching variables' = J" etc., but transformation ((7j) contains no complex- 

valued functions and is smooth over the entire spatial domain. The price for the 
smoothness is that under arctangent transformation (0), which is linear only in the 
vicinity of the origin of coordinates, familiar geometric figures become distorted on 
new mesh and, analytically, are defined by equations with x, y and z coordinates 
replaced by xotanu etc; this poses high demands on index averaging technique at index 
discontinuities for better numerical convergence, as found in Section ^ for frequency- 
domain eigenproblem. How to go around that by modifying the transformation function 
in ^ is explained in Sectional where we end up with non-absorbing PML formulation. 



4. FDFD simulation of guided modes 

It is common intuition that light tends to concentrate in high-index regions, so one may 
expect that dielectric profile transformed according to (0) shall support spurious modes 
guided along computation boundaries where some of the permittivity and permeability 
components head to infinity. To disprove this, guided propagation in step-index fiber was 
simulated with full-vector FDFD algorithm implementing permittivity and permeability 
profiles as in ((Tj), adopted to the two-dimensional geometry. No spurious boundary- 
guided solutions were found, and numeric results show reasonable agreement with those 
obtained analytically for the same fiber (figure^). Convergence was found affected by the 
scaling parameters Xq, yo (the two curves in the figure correspond to the Xq = yo = 3R 
and Xq = yo = 4i? choices, where R is the fiber radius), and very sensitive to index 
averaging scheme. This latter sensitivity is due to highly steep and nonlinear profiles of 
e and in squeezed coordinates, which can lead to systematic under- or overestimation 
of the eigenvalues if inappropriate index averaging at material interfaces (e.g., simple 
volume-weighted averaging, quite widespread in FDTD and FDFD modeling) is used. 
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Figure 1. Relative discrepancy Aricff/ricfr between SOB-FDFD computed and 
analytically calculated mode indices of step-index fiber (radius R — 2 /im, material 
index 1.45) in air background at 1.45 fim, with increasing the number of grid points 
N per domain width. 

5. A new way to construct PMLs 

An interesting question is how arbitrary the transform function in (0) is. Indeed, one 
may use, e.g., hyperbohc tangent for the mapping; such choice would lead to only slight 
changes in (jH}, 0, with cos^(M7r/2) etc. on (—1, 1) substituted by 1 — on the same 
domain. And indeed, numeric simulations in tanh-transformed coordinates show results 
similar to those in Section |3J although convergence was found poorer for the example 
in figure ^ probably owing to higher steepness of appropriately transformed e and /i 
profiles near the domain boundaries which spur numerical errors at discretization. 

Going a step further, one might introduce piecewise mapping functions like u = x for 
Ixl < Xi, u = Xi-\-X2 tanh for I a; I > Xi, where Xi defines the interface between space 
region untouched by the transformation, and X2 the width of 'non-absorbing perfectly 
matched layer'. The analogy with conventional dispersive PMLs becomes even more 
pronounced if we put X2 oc A, which is a rather natural choice as noted in Section El 
The advantage of this 'piecewise' formulation over smooth arctangent or hyperbolic 
tangent squeezing is that within the \x\ < Xi, \y\ < yi, \z\ < zi region, e and /i profiles 
are defined as on untransformed grid; the disadvantage is that they are consequently 
sharper near computation domain boundaries, for the domain of the same width. 

The notions of 'complex coordinate stretching' ^Ij or 'lossy mapping of space' ^H] 
have long been used to derive standard (lossy) PMLs, though it was not always clear 
what physical sense those complex-valued 'degrees of freedom' added by hand to spatial 
variables in Maxwell equations do have; in operational terms, whether complex-valued 
coordinates are observable. The proposed SOB technique paves the way to construct 
lossless PMLs with a clear geometric meaning of outer space squeezing; and if losses 
should be introduced (as in frequency-domain calculations of leaky modes), this can 
be done by adding, under the leakage irreversibility condition, an imaginary part to 
refractive index of surrounding medium before squeezing. 
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6. Conclusion 



Squeezing of open boundaries is proposed as an inexpensive and, at analytic stage, 
rigorous alternative to standard lossy PML technique. What makes our method so 
attractive is its conceptual clarity: we do not surround computation window with 
artificial lossy media; we do not modify Maxwell equations in any way; all we do is 
we choose coordinate system allowed by covariant nature of Maxwell equations and 
suitable for calculations — and for the finite-difference or finite-element calculations on a 
bounded domain, a suitable system is one that has bounded coordinates. The method 
is more straightforward to apply in time domain; in our proof-of-principle frequency- 
domain simulations of guided propagation, no spurious modes confined in the regions 
of strongly modified e and /i have been detected. 
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